Stochastic series expansion method with operator-loop update 
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A cluster update (the "operator-loop") is developed within the framework of a numerically exact 
quantum Monte Carlo method based on the power series expansion of exp(— /3H) (stochastic series 
expansion). The method is generally applicable to a wide class of lattice Hamiltonians for which 
the expansion is positive definite. For some important models the operator-loop algorithm is more 
efficient than loop updates previously developed for "worldline" simulations. The method is here 
tested on a two-dimensional anisotropic Heisenberg antiferromagnet in a magnetic field. 



The path integral formulation of quantum statisti- 
cal mechanics is a useful starting point for numerical 
studies of interacting many-body systems in cases where 
positive-definiteness can be assured. Monte CmwIo algo- 
rithms based on the "Trotter decomposition"Eli3 in dis- 
crete imaginary time, commonly referred to as "world- 
line" methods, have been used extensively for studies of 
quantum spins and bosons, as well as fermions in one di- 
mension (in higher dimqijsions the fermion path integral 
is not positive definite) .El Recently, two important tech- 
nical developments have lead to significantlv^ more effi- 
cient simulation algorithms. A generalizationa of ckister 
updates used in classical Monte Carlo simulationsEl can 
reduce the autocorrelation times of some simulations by 
orders of magnitude, nfl thereby enabling studies of mod- 
els in parameter regimes where standard local updating 
schemes do not efficiently explore the configuration space. 
Algorithms have also been constructed that work directly 
in the imaginary time continuum JSEj thus producing re- 
sults free of systematic errors without the extrapolations 
to zero discretization which are required in order to ob- 
tain numerically exact results using the Trotter decom- 
position. 

There are, however, still unresolved issues for these im- 
proved algorithms. For some important models the loop 
schemes do not take into account all interactions in the 
system, and hence an a posteriori acceptance probabil- 
ity has to bg [assigned after the loop-clusters have been 
constructed. Q'EJ This can seriously affect the efficiency of 
simulations. §f«?ie loop algorithms also break down due 
to "freezing" ,Q'll3 when the probability is high for a single 
cluster to encompass the whole system. It is also often 
a highly non-trivial task to construct an algorithm for a 
new Hamiltonian — it would clearly be desirable to have 
a simple recipe for an arbitrary model. 

In this Communication, a general loop-type updating 
scheme is constructed within the "stochastic series ex- 
pansion" (SSE)e[E3 framework. This approach to quan- 
tum simulations is based on sampling the diagonal ma- 
trix elements of the power series expansion of exp(— /?_ff ) 
[where H is the Hamiltonian and /3 the inverse temper- 
ature] and is related to a less general method proposed 
by IIandscomb.li3 The SSE scheme is as general in ap- 
plicability as the worldline method, and like the contin- 
uous time variant, it is numerically exact (there is also 



a strong relationship between the two methodsEiJ). SSE 
algorithms have been applied to numerous problems over 
the past several years, but so far only local updating 
schemes have been used. The "operator-loop" algorithm 
introduced here has the same favorable effects on auto- 
correlation times as the loop updates developed within 
the worldline scheme. In addition, the method overcomes 
the problems discussed above; all interactions are taken 
into account in the loop construction, there does not ap- 
pear to be any problems related to freezing, and the al- 
gorithm is very easily implemented for a wide range of 
models. 

For definiteness and sake of simplicity, the operator- 
loop algorithm will here be described for simulations of 
the anisotropic S — 1/2 Heisenberg model in a magnetic 
field, defined in standard notation by the Hamiltonian 

H = JJ2[AS:^S^ + '^{S+Sr + Srs+)] - hJ^S^, (1) 

where {i,j) denotes a pair of interacting spins on a lat- 
tice in any number of dimensions. In addition to serving 
as an illustration for a general SSE operator-loop algo- 
rithm, simulation results for this model will show explic- 
itly that problems present with other loop algorithms are 
avoided. With the standatd, worldline loop algorithms, 
freezing occurs for A > l.Qli3 The loop construction alaa 
does not take into account a non-zero magnetic field /i,E£l 
hence making simulations of large h > systems prob- 
lematic. In the present algorithm, h is explicitly taken 
into account in the loop construction and simulation re- 
sults show that A > 1 poses no problems. 

For the construction of the SSE configuration space 
the Hamiltonian is first written as 



M 



H — —J 2_^[Hi.b — H2^b], 



(2) 



b=l 



where Hi^ and H2,b are symmetric bond operators cor- 
responding to an interacting spin pair {i{b), j{b)); 



Hi,b — C 

H2,b 
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2j(^i(b) + ^i(b)) 
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2(^t{b)^3{b) 



^,.{b)^j{b)) 



(3) 



The constant C only shifts the energy and can be cho- 
sen to assure a positive definite expansion for any non- 
frustrated lattice. The number of spins in the system is 



denoted by TV; the number of bonds M = Nd for a cubic 
lattice in d dimensions. 

The partition function Z = Trje^'^-^} is expanded as 
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in the basis {|a)} = {|5f , S'f, . . . , S^)}. This expansion 
converges exponentially for n '^ N j3. A truncation at n = 
L of this order is imposed, and a unit operator -ffo,o = 1 
is introduced to rewrite Eq. fl) as (for a more thorough 
discussion, see Ref. Ill) 
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where S^ denotes a sequence of operator-indices; 
Sl = [ai,bi]i, [a2,62]2, • ■ • , [olj^lIl, 



(5) 



(6) 



with Oi G {1,2} and h G {0,...,M}, or [ai,6i] = [0,0], 
and n denotes the number of non-[0,0] elements in S^. 
In principle, each term in (||) should be multiplied by 
a factor (—1)"^, where n2 is the total number of [2,6] 
elements in Sl- However, for a non- frustrated lattice 
this number must always be even for the matrix element 
to be non-zero. Choosing C in (H) such that all matrix 
elements of -ffi.f, are positive, the expansion is then pos- 
itive definite. A Monte Carlo procedure can therefore 
be used to sample the terms. £a, S*!,) according to their 
relative weights. Previouslyjj'til sampling schemes were 
devised based on (i) local substitutions of single diago- 
nal operators, [0,0]p <-> [l,6]p, and (ii) pairs of diagonal 
and off-diagonal operators [1, h\p^ [1, b]p^ <-> [2, b]p^ [2, b]p^. 
The diagonal update (i) will also be used here. The 
new operator-loop update involves any number of diago- 
nal and off-diagonal operators and is much more efficient 
than the simple pair update (ii). 

It is convenient to introduce the notation \oi{p)) for 
states obtained by acting on \a) in Eq. (|^) with the first 
p elements of the operator string. 



|a(p))-[]i/a,,6.|a) 



(7) 
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and to define states \ab{p)) = l'S'f(f,)(p), 5'J(jj(p)) 
bonds. For a contributing term, \a{L)) = |a(0)) — \a). 

The simulation starts with some random state |a) and 
an operator string [0, 0]i, . . . , [0, 0]l containing only unit 
operators. The cut-off L is chosen arbitrarily and ad- 
justed during the equilibration phase of the simulation so 
that it will always be larger than the highest n reached 
(hence leading to no detectable truncation error) . The di- 
agonal update [0,0]p <-> [l,6]p is carried out sequentially 
at each position p — 1, . . . , L for which [op, bp] — [0, 0] 
or [1,6]. When accepted, such an update changes the 
expansion power n by ±1. Acceptance probabilities that 



satisfy detailed balance are obtained using Eq. (^) and 
the fact tiiat there are M random choices for 6 in the -^ 
direction;llJ 



P([0,0]p->[l,6]p) 
P([l,6]p^[0,0]p) 



M(3{ab{p)\H^,b\ab(p)) 

L — n 

L-n + 1 

MrJ{ab{p)\Hi,b\ab{p))'' 



(8) 



where a number larger than 1 should be interpreted as 
probability one. The state |a(0)) is stored at the be- 
ginning of an updating cycle. Each time an off-diagonal 
operator [2, b]p is encountered, the corresponding spins 
are flipped so that the states in Eqs. (||) will be available 
when needed. 

The second, new type of update is carried out with 
n fixed. It is then convenient to disregard the [0, 0] unit 
operator elements in Sl and instead work with sequences 
Sn containing only the Hamiltonian operators [1,6] and 
[2, 6]. The propagation index p will in the following refer 
to this reduced sequence. Further, full bond operators 
including both the diagonal and off-diagonal terms are 
defined; Hi, =^ Hij, + H-z.b- The matrix element in dq) 
can then be written as 

n 

M(a,5„) = l[{ab^ip)\HbJab^{p-l)). (9) 

p=i 

The non-zero matrix elements are 

{iA\Hb\i,i)=C-A/A-h/{2J), 

(T,TimiT,T) = C-A/4 + V(2J), 

(i, T l^;,! i, T) = (T, i \Hb\ ti) = C + A/4, (10) 

(T,i|i?;>li,T) = (i,TI^6|T,i)-i/2. 

C should be chosen such that all diagonal matrix ele- 
ments are larger than (or equal to) zero. M(a, Sn) can 
be graphically represented as a set of n vertices connected 
to the propagated spins, as shown in Figure 1(a) for a sys- 
tem with 4 spins. Two spin states "enter" each vertex, 
and "exit" in either the same states or flipped (the direc- 
tion of the propagation is clearly irrelevant). The allowed 
types of vertices, corresponding to the non-zero matrix 
elements ([lO|), are shown in Fig. 1(b). Fig. 1(a) displays 
all the full spin states at each "event" , but clearly there 
is much redundant information in this picture. The spin 
states at the four "legs" of the n vertices completely spec- 
ify the full spin conflguration (except for spins that hap- 
pen not to be connected to any vertex). In order to carry 
out the operator-loop update, a linked list of the vertices 
with their four spin states is constructed using the cur- 
rent state \a) and the index sequence Sl- The list is 
doubly linked, so that it is possible to move in either di- 
rection from any leg of a given vertex to the leg of the 
next or previous vertex connected to the same spin. 

The principles of the operator-loop update are now 
quite simple to state: One of the n vertices is flrst chosen 
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FIG. 1. (a) Representation of a matrix element product 
M{a,S„), Eq. (H), with n = 7, for a 4-spin system. The 
vertical solid and dashed lines indicate the spin states acted on 
by the operators Ht, which are represented by the horizontal 
bars, (b) shows the allowed vertices, which are associated 
with the non-zero matrix elements (urn. 
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FIG. 2. The four paths on a vertex in the case of the en- 
trance point being the low-left leg. The entrance and exit legs 
are indicated by the arrows. The spins at these legs are flipped 
in the process; the states at the other legs remain unchanged. 



at random, and one of its four legs is randomly selected 
as the entrance point. One of its legs is then chosen as 
the exit point from the vertex, according to probabihties 
to be specified below. The four possible vertex paths, 
in the case of the entrance being the low-left leg, are 
illustrated in Figure 2. The spins at both the entrance 
and exit legs are flipped. [Note that the entrance and 
the exit can be the same leg. Fig. 2(d), in which case the 
net effect is no spin flip; only a reversal of direction of 
movement in the list.] The chosen exit leg points to a 
leg of another vertex in the linked list, the spin at which 
is also flipped. From this vertex, an exit leg is again 
chosen, which points to another vertex, e.t.c. After some 
varying number of steps, the exit of the last visited vertex 
will point to the original entrance point of the update. 
The loop then closes and the result has been to flip all 
the spins along the random path followed in the process. 
Since the operator list is a periodic structure (because 
\a{n)) — |a(0))), any state \a{p)) can be affected in the 
update, and the sum over states \a) in Eq. (pi) is therefore, 
implicitly, also sampled in the process. 

The probabilities for the four different choices of exits 
from a given visited vertex are simply proportional to 
the matrix elements (M) corresponding to the resulting 
vertices, i.e., those where the spins at the entrance and 
exit legs have been flipped. It is intuitively clear that this 
operator-loop procedure satisfies detailed balance and, in 
combination with the diagonal single-operator update, is 
ergodic in the grand canonical ensemble (fluctuating total 
z-component of the magnetization) , including all winding 
number sectors. For lack of space, a rigorous proof will 
not be presented here. 



Note that one of the paths (a)-(c) in Fig. 2 will al- 
ways have zero probability, since the Hamiltonian (|l|) 
does not contain operators S^ S^ or S~ S^ . These op- 
erators could be included in a more general model and 
then all four paths would be allowed. The probability of 
the "bounce" process (d) is always in principle non-zero. 
However, in some cases it is possible to exclude this path. 
Consider the XY model in zero field, i.e., A = h = 0. If 
C = 1/2 is chosen, all the non-zero matrix elements in 
( pl]| ) equal 1/2. Detailed balance is then satisfied also by 
only choosing, with equal probabilities, among the two al- 
lowed paths (a)-(c). For the isotropic Heisenberg model, 
i.e., A = 1, h = 0, and with C = 1/4, the bounce can also 
be neglected. The only allowed path is then always the 
"switch and reverse" (c) [which corresponds to a substi- 
tution [1, b] ^^ [2, b] in terms of the operators in Sl], and 
hence the loop construction is completely deterministic 
in this important case. 

A full updating cycle consists of the following steps: 
First the diagonal single-operator update is carried out 
at all positions in Sl with diagonal operators. The linked 
list of vertices is then constructed and a number of loop 
updates are performed. The typical size of a loop de- 
pends strongly on the model parameters. The number of 
loops to be constructed in each cycle is therefore chosen 
such that on average a total of ~ (n) vertices are vis- 
ited. The updated vertices are finally mapped onto the 
corresponding operator- indices [a, b] and written into Sl- 

To demonstrate the efficiency of the new algorithm, re- 
sults are next presented for two different cases where, pre-, 
vious loop algorithms have encountered difficultiesO'E^ 
The anisotropic model in zero field and the isotropic case 
with a field. The estimators for various observables of in- 
terest have been discussed in detail in Ref. |lj. The cor- 
rectness of the simulation code was verified by comparing 
results for a 4 x 4 lattice with exact results obtained by 
diagonalizing the Hamiltonian. The results to be pre- 
sented next were obtained using lattices sufficiently large 
to eliminate finite-size effects. For the lowest tempera- 
tures considered, 64 x 64 spins were typically used, and 
on the order of 2 x 10^ updating cycles were carried out. 

The susceptibihty, x = f^jSEz S-f)/^^ for the case 
h — \s shown in Figure y for several values of the 
anisotropy_A, Unlike with the standard worldline loop 
algorithm,Qll3 there are no problems with "freezing" in 
simulations for A > 1. The exponential decay of x to 
as T — > for A > 1 reflects the opening of a gap 
in the spectrum for these systems. For the isotropic 
case (A — 1), the results, are in perfect agreement 
with previous calculations Il3 For the XY model (A = 
0), a temperature-independent behavior is seen at low- 
temperature {T / J ^ 0.2), in agreement with a predic- 
tion of chiral perturbation theory.llZl Quantitatively, the 
T- independent value should be x = Ps/c^,lII where ps is 
the spin stiffness and c the spin- wave velocity. The result 
X — 0.2095(3) obtained here at T/ J = 0.05 is consistent 
with this priediction and recent ground state calculations 
of Ps and c.ta 
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FIG. 3. Magnetic susceptibility vs temperature for the 
zero- field Heisenberg model with anisotropy parameter A = 0, 
0.5, 1, 2, and 4 (top to bottom). 
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FIG. 4. Magnetization vs temperature for the isotropic 
Heisenberg model (A = 1) in magnetic fields h/J = 2, 1, 
1/2, 1/4, and 1/8 (top to bottom). 



The magnetization per spin, m = (^^ Sf)/N, is shown 
for an isotropic interaction and several strengths of the 
magnetic field in Figure ^. For all field-strengths, there 
is a maximum in m between T/J = 0.5 and 1, reflect- 
ing the cross-over between high-temperature independent 
spin behavior and antiferromagnetic correlations devel- 
oping at lower T (also seen in the zero-field susceptibility 
in Fig. ph. Note the shallow minimum at lower temper- 
atures for h/J < 1. This reflects the temperature scale 
at which the local, short-range antiferromagnetic corre- 
lations are the strongest. 

The operator- loop simulations are very efficient for any 
strength of the field, since a /i > is taken into account 
in the loop construction. With other loop algorithms ,&t£l 
an a posteriori acceptance probability has to be assigned 
for updates in which the total magnetization changes. 
This acceptance probability decreases rapidly with in- 
creasing field strength, leading to an autocaurelation 
time which increases exponentially with h/T.tJ Previ- 
ous simulationst3 were therefore restricted to h/T < 4. 
Fig. Q shows results up to h/T — 40, and there are no 
signs of increasing autocorrelation times even for much 
higher values. 

To conclude, the operator-loop algorithm introduced 
here has several-advantages over other loop methods sug- 
gested recently.crEJ The most important is that all inter- 
actions, including external fields, are taken into account 
in the loop construction, thus eliminating the need for a 
posteriori acceptance probabilities_that restrict the appli- 
cability of the previous methods.Ej Like the [Cpntinuous- 
time version of the worldline algorithm,El"U the SSE 
method is completely approximation free. The config- 
uration space is, however, discrete, and the only floating 
point operation required in the simulation is the gener- 
ation of uniformly distributed randoni^iiumbers. In the 
continuous-time worldline algorithms,crllll on the other 



hand, high-precision values of imaginary times have to 
be manipulated. One can therefore expect that the 
operator-loop algorithm is faster in many cases, in par- 
ticular for the uniform Heisenberg model, where the loop 
construction is deterministic. It is also interesting to note 
that certain expectation values have simpler cstiroators 
in the SSE framework than for worldline methods .Eil 

The method has here only been demonstrated for the 
anisotropic Heisenberg in a magnetic field. Generaliza- 
tions to other models with two-body interactions are al- 
most trivial, however. The vertices depicted in Fig. 1 
only involve other degrees of freedom at the "legs" . For 
example, for Hubbard-type models the legs can have 
charge c = 1 and spin s = ±2,ors = and c = 0, 2. The 
vertex paths in Fig. 2 then involve changing these quan- 
tum numbers by some values {Sc, 6s) at the entrance leg, 
and changing them by (Sc, Ss) at an exit leg in the same 
direction [paths (a) and (b) in Fig. 2] or {—Sc, —Ss) at 
an exit in the reverse direction [paths (c) and (d)]. Im- 
plementation for a new model thus essentially involves 
specifying all allowed vertices, i.e., all non-zero matrix 
elements of type (|l^) . 
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